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Abstract. The following electromagnetism (EM) inverse problem is addressed. It consists 
in estimating local radioelectric properties of materials recovering an object from the global 
EM scattering measurement, at various incidences and wave frequencies. This large scale ill- 
posed inverse problem is explored by an intensive exploitation of an efficient 2D Maxwell solver, 
distributed on High Performance Computing (HPC) machines. Applied to a large training data 
set, a statistical analysis reduces the problem to a simpler probabilistic metamodel, on which 
Bayesian inference can be performed. Considering the radioelectric properties as a dynamic 
stochastic process, evolving in function of the frequency, it is shown how advanced Markov 
Chain Monte Carlo methods, called Sequential Monte Carlo (SMC) or interacting particles, can 
provide estimations of the EM properties of each material, and their associated uncertainties. 



1. Introduction 

The inverse problem is described in figure [TJ The Radar Cross Section quantifies the scattering 
power of an object, at a given incidence and wave frequency. It is defined as the ratio 
between the radar transmitted power and the incident power density (in plane wave) [Tj. RCS 
measurement process is schematically presented on the right part of figure [1} The object or 
mock-up is illuminated by a quasi-planar monochromatic wave, inside an anechoic chamber 
where interferences are limited. The acquisitions are realized at K successive discrete frequencies 
(/i)/2) - "" >/k)> f° r different incidence angles (by piloting motorized rotating support). From 
this raw data, a signal processing is performed, mainly consisting of calibration and filtering. 
At the end, the measurement provides an evaluation of the calibrated complex (amplitude and 
phase) scattering coefficient, for each frequency and incidence. 
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Figure 1. The EM inverse problem 



A metallic axi-symmetric object is recovered with TV" areas (see figure [3]), each area 
corresponding to a material with its associated isotropic radioelectric properties, i.e. the complex 
parameters of permittivity e and permeability /x. The EM inverse problem can be expressed as: 
is it possible to extract some local information on the material properties (e^,^) of each area 
from the global scattering measurement? 

2. The stochastic modeling 

At a given frequency fk, the system state X& can be defined by (omitting fk to lighten the 
notations): = [e ; e" fi' [x"] T ', with e', e", /j,' and /i", respectively the real and imaginary 
permittivity and permeability components of the N areas. Notice that the state dimension can be 
high (4- N with N ~ 100). On the other hand, considering measurements at the given frequency 
fk with various incidence angles in both polarizations (TM and TE), the observation vector is 
made of the real (3?(-)) an d imaginary (3 : (-)) parts of the complex scattering coefficients Cxm 
and cte measured at M angles 9i,--- ,9m- Y& = [3J(ctm) 9(cte) S(ctm) 9(cte)] T - 
Considering all frequencies, vectors X = (Xi, . . . , X«-) and Y = (Yi, . . . , Y#) respectively 
define complete system state and observation. 

2.1. Observation model 

The 2D-axisymetric Maxwell solver software (J~2d) can predict the observation from the system 
state. Assuming a multidimensional Gaussian measurement uncertainty model, it leads to the 
following likelihood model Y& | X& ~ A/"(J r 2£)(Xfc), Rfe) (at a given frequency fk), where Rfc is 
the covariance matrix (assumed known). To avoid numerous and heavy T%b computations, we 
have developed the following global approach. First, the high dimension state space and the 
associated system response are explored randomly around expected properties (prior knowledge), 
computations being massively distributed on HPC machines. 

Let B A = {(Xj^, y£ 1} ), • • • , (X.[ Ne) , Y ( k NE) )} the training data composed of N E couples. It is 
then processed by N-D statistical techniques; sensibility analysis and model reduction techniques 
can possibly reduce the state space dimension. Applying multidimensional regression, it turns 
out that the model is approximatively linear. According to the studied cases, the linearity errors, 
evaluated by residue analysis and bootstrap techniques, are significant but much lesser than the 
RCS measurement uncertainties. Finally, the following linear Gaussian model can be considered 
(at a given frequency fk)- 

Y fc |X fc ~AA(A fc -X fe +Y°,R fc ) (1) 

where the deterministic part of the linear model is given by the learned matrix 4M x 4A A^ 
and the vector Y°. 

2.2. Prior model at a fixed frequency 

Let us fix a frequency fk- We model our a priori knowledge on Xj, with a Gaussian distribution 
A/"(rrifc,Pfc). The object is divided in iV& blocks of areas, each of them composed of a rather 
homogeneous material. The location of these blocks is known exactly. Prior mean value 
is defined with reference values of e', e" , // and fj," for each of theses blocks. Then, for any 
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component XI, we define a variance o~ l Sk as a mix between absolute and relative uncertainty. 
To take into account the spatial local homogeneity, covariance matrix is defined block by 
block independently (and separately for each of the e', e" , fjf, fx") by correlation relations between 
Xfc's block-sharing components: 

Cav(KlX{) = at/ ■ 



where ps E [0, 1] is a spatial correlation parameter (typically ps = 0.95 in our applications). 
It means that the correlation between 2 areas dicreases geometrically at speed ps with the 
distance between them. 

2. 3. Inter- frequential prior model 

Radioelectric properties e', e", p! ', p" are known to vary in function of the wave frequency [lj; 
their non-stationary dynamic can be quite different from frequency /i to fx- However, in order 
to take account of expected frequency profiles regularity (for each material and e', e", p' , p"), 
we model sequence (X&, k G {1, . . . , K} with a linear Gaussian correlation structure, given by 
the generalised autoregressive (AR) random process [2]: 

Xi ~ AA( mi ,Pi) 

(X fe+1 - m fc+ i) = D p • H fc+1 • H^ 1 • (X fc - m fe ) + y/l d - ■ H fc+1 • V k (2) 

where 

• for each k, H/% is the "square root" matrix of covariance matrix Pj., defined as being the 
unique symetric definite positive matrix satisfying • = P&; 

• (V fc , k G {1,...,K}) are i.i.d. Af(0, I d ); 

• Dp is a positive diagonal matrix commuting with the H&, further described, and depending 
on a parameter p. 

One notices that so defined sequence X^ still admits M (m&, P^) as marginal distributions. More 
generally, it can be shown that the distribution of concatenated vector X = (Xi, . . . , X#) is a 
Gaussian distribution with mean m = (mi, . . . , mjf) and covariance matrix: 
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where T~L is the diagonal by block matrix T~L = diag(Hi, . . . , H^-), which explicits every joint 
distribution (Xi, Xj). 

Finally, let us clarify parameter p (and matrix Dp). D p 's role in these equations makes clear 
that it's a frequential correlation parameter. Its dimension is to be chosen among 3 possibilities 
according to our assumptions. 

• 1 st case : p is assumed independant from material and e', e", p', p" . It is one-dimensional, 
p £ [0, 1], and D p = p.I d . 

• 2 nd case : p depends on the material (block). It is iV^-dimensional, p G [0,1]^, and Dp 
giving to each line (component of X^) its associated component of p. 

• 3 rd case : p depends on material and e', e", p', p". It is 4. ^-dimensional, and D p is the 
diagonal matrix composed with 4.A^ pj.I d -type blocks. 



2-4- Global model 

Conditionnaly to frequential correlation parameter p, the problem of the determination of 
X = (Xi, • • • , X#) given the measurements Y = (Yi, • • • ,Yjf) can be expressed as a classic 
linear Gaussian hidden dynamic Markov process observed at "times" fk (k = 1, • • • , K): 

X fe+1 = M£ • X fe + and Y fe = A k • X fc + wf 

where wjP""' and w[ are independant Gaussian noises with known parameters, and a known 
matrix depending on p (see ([I]) and Q). Parameter p, intuitively representing the inter- 
frequency regularity, is unknown and to be estimated. In respect to the probabilistic point 
of vue, it is probabilized, and given a prior distribution p{p). 



3. Sequential Monte Carlo approach for global inversion 

3.1. Rao-Blackwellised SMC algorithm 

The posterior distribution p(X, p|Y) can be decomposed as: p(X, p\Y) = p(X.\p, Y) -p(p\Y). 
Conditionally to p, the system is linear Gaussian: the conditional distributions p(X|p, Y) can 
be straightforwardly computed by Kalman filtering, including in this off-line context backward 
Kalman smoothing. On the other hand, the term p(p\Y) oc p(Y\p) ■ p{p) can be evaluated 
(up to a normalising constant) for a given p using the likelihood term provided by the Kalman 
filter and prior distribution p(p). Consequently, in order to exploit this conditional structure of 
the system, Kalman smoothers are applied and integrated in an interacting particle approach. 
This idea of mixing analytic integration (here Kalman evaluation of p(X.\p, Y)) with stochastic 
sampling is a variance reduction approach, known as Rao-Blackwellisation [3]. 

Similarly to [3], we choose to implement an efficient interacting particle approach, in order 
to estimate the marginal distribution r/(p) := p(p\Y) oc p(Y\p) ■ p(p). Sequential Monte 
Carlo (SMC) is a stochastic algorithm to sample from complex high-dimensional probability 
distributions. The principle (e.g., |4j) is to approximate a sequence of target probability 
distributions (i] n ) by a large cloud of random samples termed particles ((n)i<k<N p £ E Np , 
E being called the state space. Between "times" n — 1 and n, the particles evolve in state space 
E according to 2 steps: 

• a selection step: every particle C n -\ ls given a weight uji defined by a selection function g n 
(uji = g n {C n _i)). By resampling (stochastic or deterministic), low- weighted particles vanish 
and are replaced by replicas of high- weighted ones. 

• a mutation step: each selected particle C n -\ m ove, independently from the others, 
according to a Markov kernel M n . 
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Evolving this way, the cloud of particles, and more precisely the occupation distribution 
Vn p = ~^~^2k=i^Ck ( sum °f Dirac distributions), approximates for each n the theoretical 
distribution r\ n defined recursively by the Feynman-Kac formulae, associated with the potentials 
g n and kernels M n [5]. 

Back to our objective of sampling from rj, we then define the sequence of distributions f] n (p)'- 



r] n (p) (xp(Y\p) an -p(p) 

where {a n )o<n<n f is a sequence of number increasing from to 1, so that: rjo is prior 
distribution p(p), easy to sample, r\ nj is target distribution rj and sequence (<r) n ) admits a 
Feynman-Kac type structure with calculable selection functions g n and Markov kernels chosen 
so that r\ n .M n = M n (Metropolis-Hastings for example). The distribution rj is then interpreted 
as being the last distribution of a Feynman-Kac sequence, on which SMC can be performed, the 
estimator of rj being the occupation distribution r] Np extracted from the last cloud of particles. 

3.2. Results 

The occupation distribution r\ p which approximates r] = p(p\ Y) can be represented dimension 
by dimension via histograms (see figure [2J . 




Figure 2. Estimation of a 5-dimensional p's distribution with N p = 100 particles. 

For each frequency f k , this approximation rj Np ~ r], associated with the theoretical 
conditioning relations 

E(X k \Y) = E[E(X k \p,Y)\Y] 

Cov(X k \Y) = E [Coa(X k \p, Y)|Y] + Cov [E(X k \p, Y)|Y] 

can deliver estimators of respectively the mean and the covariance matrix of .Roughly 
speaking, the posterior estimation is performed by randomly picking a pi from the final cloud 
of particles and computing associated samples of by a Kalman smoother conditionally to p{. 
It is illustrated in figure [3j with a good agreement between the true state and estimated state. 

Moreover, for any fixed area, the method provides estimators of the mean and marginal 
variance for every frequency, so that the results can be presented as frequential profiles, with 
marginal uncertainties (see figure [4]). Even when the true (simulated) EM property profiles are 
chosen markedly divergent from the prior AR-type model, the method turns out to be robust. It 
results from the adaptive estimation of p which provides small values of p (i.e. weak correlation 
of EM properties for close frequencies) in the case of highly irregular true profiles. 

4. Conclusion 

An efficient statistical inference approach has been applied to estimate local material 
radioelectric properties from global EM scattering measurements. It combines intensive 
computations, meta-modeling and advanced sequential Monte Carlo techniques dedicated to 
frequency dynamic estimation. 
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Figure 3. Estimated radioelectric properties of a Nf, = 5-block object (N = 19 areas), at fixed 
frequency f k . 
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Figure 4. Estimated radioelectric properties (K = 30 frequencies) of a fixed area. 
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